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Abstract 

Granular materials will segregate by particle size when subjected to shear, as occurs, for 
example, in avalanches. The evolution of a bidisperse mixture of particles can be modeled 
by a nonlinear first order partial differential equation, provided the shear (or velocity) is a 
known function of position. While avalanche-driven shear is approximately uniform in depth, 
boundary-driven shear typically creates a shear band with a nonlinear velocity profile. In this 
paper, we measure a velocity profile from experimental data and solve initial value problems that 
mimic the segregation observed in the experiment, thereby verifying the value of the continuum 
model. To simplify the analysis, we consider only one-dimensional configurations, in which a 
layer of small particles is placed above a layer of large particles within an annular shear cell and 
is sheared for arbitrarily long times. We fit the measured velocity profile to both an exponential 
function of depth and a piecewise linear function which separates the shear band from the rest 
of the material. Each solution of the initial value problem is non-standard, involving curved 
characteristics in the exponential case, and a material interface with a jump in characteristic 
speed in the piecewise linear case. 

1 Introduction 

When set in motion through vibration or shear, granular materials have a strong tendency to 
segregate into bands containing particles of similar size, shape, or density [12] . In this paper, 
we focus on shear-induced segregation by size, which appears in a variety of configurations and 
applications, including avalanches [T3], rotating tumblers |lUj . and internal shear experiments [S]. 

Continuum models of avalanche flow have been derived using ideas from shallow water theory, 
in which a thin-layer approximation captures both the free surface shape and the underlying depth- 
averaged velocity |15j . These models typically do not account for segregation. However, segregation 
in avalanching flows has been modeled by a mass transport equation alone, in which a roughly 
constant shear rate is specified [U [16]. Here, we adapt a mass transport segregation model to 
situations where the shear rate is far from constant, reflecting the nonlinear dependence of the 
velocity of particles on position, as is commonly the case for boundary-driven flows jTT]. Through 
experiments on a mixture of two particle sizes within an annular shear cell, we measure the velocity 



1 



window 



Figure 1: Initial configuration of particle layers in the experimental annular Couette cell. 



profile and incorporate the resulting spatially-dependent shear rate into the constitutive law of the 
model. 

The Gray-Thornton model jlj for segregation by size in an avalanche containing two species of 
particles with similar density but different size, takes the form 

tpt + u(z)(p x + (w(tp,z)(p) z = 0. (1.1) 

In this partial differential equation (PDE), ip(x, z,t) represents the concentration (fraction by vol- 
ume) of smaller particles as a function of the distance down the avalanche x, distance z above the 
base and time t. The bulk flow is represented by the velocity u(z) parallel to the base; the normal 
velocity w((p,z) of small particles is due to segregation dynamics. Both u(z) and w(ip,z) are as- 
sumed to be known functions whose functional forms are deduced as part of the model derivation. 
In avalanche flow, u(z) is roughly linear near the surface [11], so that the shear rate |tt'(.z)| is close 
to constant, and w(cp, z) = — k(l — cp) for positive, constant k proportional to the constant shear 
rate. Thus, in the Gray-Thornton model, the segregation rate is independent of depth. Note that 
1 — 99 is the concentration of large particles, so that this form for w(<p,z) may be considered to 
represent the availability of large voids created by relative motion of large particles. 

In this paper, we are interested in the influence of non-uniform shear rate |u'(z)|, for which 
the segregation rate will be different at different depths. For example, if there is no shear, then 
there should be no tendency towards segregation, whereas a large shear rate should induce rapid 



segregation. Our model fits into the general framework of equation (1.1), but the normal speed 
w of small particles depends on z through the depth-dependence of the shear rate |«'(2r)|. In our 
model, we assume that w is proportional to |?/(<z)|. In principle, w could be any increasing function 
of shear rate. 

To simplify matters in both the experiment and the model, we begin with a bidisperse granular 
material in which the two sizes of particles (with the same density) are arranged in a one-dimensional 
configuration, as shown in Fig. [I] It is reasonable to assume that in the subsequent evolution from 
this normally-graded configuration to an inverse-graded configuration, the concentration of each 
size of particle at each location depends only on depth and time. We explore two PDE models, 
both motivated by the structure of the velocity profile taken from experimental observations and 
differing only in the choice of depth-dependent shear rates chosen to approximate the experimental 
results. 
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The general form of the PDE we consider is the conservation law 

cpt + {sa(z)f(ip)) z = 0, 0<z<l, t>0. (1.2) 

The convex flux / : [0, 1] — > K satisfies /(0) = /(l) = 0, corresponding to the physical limits of the 
dependent variable <p. The nonconstant coefficient a{z) is the shear rate |u'(z)|; it is a monotonically 
decreasing function of the position variable z. The segregation rate parameter s > sets the time 
scale for the evolution of ip. 

We will be concerned with initial boundary value problems, with initial data corresponding to 
the one-dimensional experimental configuration: 

{0, < z < zo, 
(1.3) 
1, z < z < 1, 

and boundary conditions 

¥>(0,t) = l, <p(l,t) = 0. (1.4) 



In the experiment, shown schematically in Fig. [TJ we place a layer of small glass spheres over 
a layer of larger spheres within the annular region between fixed rigid concentric cylinders. The 
aggregate is sheared by rotating the lower confining plate at fixed vertical position. An upper 
heavy confining plate is allowed to move vertically to accommodate changes in volume, due to 
both Reynolds dilatancy |14j and changes in packing density arising from the mixing/segregation 
process [3]. The particles initially mix and then re-segregate through a process known as kinetic 
sieving: as the shearing proceeds, large particles roll and slide over one another, opening up gaps 
for the smaller particles to fall into. The small particles also act as levers for the large particles, 
which consequently tend to move vertically upwards, a process sometimes called squeeze expulsion. 
Kinetic sieving was modeled in avalanche flow by Savage and Lun [16], and subsequently by Gray 
and Thornton, using a different approach [3]. In these models (which are valid in several space 
dimensions), the segregation rate is assumed to depend only on the concentration of small particles. 
This approximation is suitable for free-surface avalanches, where shearing is provided by the effect 
of gravity, a body force. However, in our experiments, shearing is instead provided by motion at 
the lower boundary, and this is transmitted through the granular material only by particle-particle 
contacts. The resulting shear rate drops off dramatically within a few layers of particles. 

We model the depth-dependence of the shear rate a(z) in two ways, suggested by the experi- 
mental data: 

Case I: Piecewise constant shear rate: 

{k , < z < z c , 
(1.5) 
fci, z c < z < 1, 

with ko > ki > 0. 
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Case II: Smooth shear rate: 



a( z ) = a e~ z/x , < z < 1. (1.6) 

Case I is based on the observation that there is a higher shear rate near the bottom plate, 
reflecting localization within a shear band. Modeling this higher rate as a constant is a coarse 
approximation to the experimental data. However, the split into two regimes, with a material 



interface at z = z c is justified by the data. Equation (1.2) with a discontinuous function a(z) does 
not fit into the existence theory of Kruzkov [8] , and indeed, the issues of existence and uniqueness 
for this type of equation have been addressed in some generality only recently [2]. In this case, 
characteristics are straight lines, on which ip is constant, but both the characteristic speed and (p 
experience a jump at z = z c . 

For smooth functions a(z) (Case II), the existence result of Kruzkov [8] for initial value problems 
can be adapted to the initial boundary value problem, by extending a(z) and the initial conditions 
beyond the boundary: a(z) = a(0),(p(z,0) = 1, z < 0; a(z) = a(l),(p(z, 0) = 0, z > 1. However, 
characteristics are curved, and moreover, ip is not constant on characteristics. Consequently, al- 
though the structure of solutions can be characterized, the solutions cannot be found explicitly. 
The exponential form in Case II is consistent with other studies of sheared granular materials |llj . 
and provides a close fit to our experimental data. We begin the analysis of Case II by considering 
general functions a(z) that are smooth, positive and decreasing, but it turns out that the choice of 
the exponential form is particularly useful for calculating explicit solutions. 



Since our objective is to mimic the experiment, we restrict attention to initial conditions (1.3) 
that reflect experimental conditions and for which we can analyze the solutions. We compute 
these solutions in detail using the structure of hyperbolic waves. Quantitative comparison of the 
theoretical solutions of this paper with experimentally-observed segregation is presented in [9]. 

The constants z c < zq and fco > k\ > in Case I, and positive constants ao,A in Case II, 
are determined from an experimentally-measured velocity profile, and an overall segregation rate 
constant s sets the time scale. Solutions based on the experimentally determined constants are 
shown in Fig. [2j in which s is chosen to make the final time t* = 1 in Case II. The solutions involve 
a rarefaction wave, centered at (z,t) = (zq,0), in which <p(z,t) varies continuously between tp = 
and (p = 1. As the leading edge reaches z = 1 (at time t = ti), a shock wave is reflected, with a 
layer of large particles (ip = 0) growing behind it. Similarly, as the trailing edge hits the boundary 
z = (at time t = to), a layer of small particles (cp = 1) develops behind the reflected shock. The 
two shocks eventually meet at time t* , at which time the solution becomes a stable stationary shock 
representing a layer of large particles above a layer of small particles separated at z = z* = 1 — zq. 

In £j2j we describe the annular shear cell experiment and explain how we determine the model 
parameters for the two chosen cases (1.5), (1.6). In fj3| we construct solutions in each of the two 



cases. Interestingly, in Case I, the time to full segregation is independent of the shear rate ko within 
the shear band. We conclude with a discussion in £j4| 



2 Experimental Results 

The experimental configuration is an annular Couette cell (see Fig. [TJ with channel width 3.8 cm 
bounded by concentric aluminum cylinders with inner and outer radii 25.5 cm and 29.3 cm, re- 
spectively. The rotating bottom plate and an upper confining plate each have rubberized surfaces 
to enhance friction with the particles. A motor drives the bottom plate at a constant rotation 
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rate of approximately 3 revolutions per minute. The cell is filled with a 2 kg layer of glass spheres 
(diameter 3 mm) , placed above a 2 kg layer of larger glass spheres (diameter 6 mm) . The fill height 
is approximately 4.1 cm, and a heavy top plate confines the particles but is free to move verti- 
cally to accommodate changes in volume as the aggregate dilates, mixes and segregates. Further 
experimental details are available in 0|9]. 

The apparatus has a window in the outer wall, permitting us to track particle positions over time 
with a high speed (450 Hz) digital camera. The camera collects digital images at discrete intervals 
throughout the duration of the experiment, allowing us to compare particle velocities at different 
stages of the experiment. In each image, we locate the center of each particle, distinguishing 
large from small, and record the positions of individual particles. Through an automated process, 
we identify the same particle in successive frames, generating a list of the horizontal and vertical 
coordinates of each particle at a sequence of times. We refer to this list as a single-particle trajectory. 

For each single-particle trajectory, we calculate the instantaneous horizontal velocity of the 
particle as follows. First, the vertical dimension of the sample is divided into twenty-three bins 
centered at positions z = z%, i = 1, ... ,23. Each trajectory is assigned to a bin Zi based on the 
average vertical position of the particle. Using a moving interval of duration At, we determine the 
instantaneous velocity by fitting a linear function to the horizontal coordinates within At. For each 
bin, we choose an appropriate, speed-dependent, value for At, which varies from approximately 
0.1 seconds near the bottom plate to approximately 0.4 seconds near the top plate. This process 
yields a range of velocities observed for an ensemble of different particles at different times. We 
fit a parabola to the peak of the probability distribution within each bin to calculate a velocity Ui 
representing the horizontal speed of particles in bin i. 

Fig.J^a) shows the velocities m, i = 1, . . . , 23, which we refer to as the measured velocity profile. 
In the figure, we have normalized z to [0, 1] over the region of interest (described in { 2.1 ), and scaled 
the velocity so that u = 1 at z = 0. The error bar through each point (ui, zi) is the width of the 
parabola at a height one half of the maximum height, to give a sense of the distribution of observed 
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(a) 



(b) 



Figure 3: Measured velocity profile Ui (•) for (a) full cell height, with boundary layers above and 
below the dashed horizontal lines and |(b)| within the region z = [0, 1]. The dotted line is the fit to 
Case I; the solid line is the fit to Case II, as described in |2.1| The velocity profile is scaled so that 
u(0) = 1. 




(a) (b) 

Figure 4: Shear rate \du/dz\ (•) for |(a)| full cell height, with boundary layers above and below the 
dashed horizontal lines and |(b)| within the region z = [0,1]. Vertical dotted lines are the fit to 
Case I; solid line is the fit to Case II, as described in §2.1| 
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values in each bin. In total, the measured velocity profile is based on processing particle positions 
from approximately 7 x 10 5 images. In {2.1, we use the measured u% to generate appropriate 



parameters for the shear rate a(z) for use in the model. Further details concerning the collection 
and processing of the experimental data are described in [9] . 

During the processing of the data to generate the measured velocity profile, we established two 
properties which are crucial for the continuum model: 

• Velocities are similar for both large and small particles; calculating m separately for large 
and small particles gives negligible differences. 

• Velocities reach steady-state after a short initial transient of approximately 0.05 t* xp , where 
i* xp = 700 seconds is the duration of the experiment. This observation justifies the use of a 
time-independent velocity profile u(z) in the model. 

2.1 Determining the Shear Rate Profiles 

To obtain the position-dependence of the shear rate from the measured velocity profile, we first 
take finite differences of u\ between adjacent layers z%. 

u'{z i)xi Ul ± i - — — , Azi = Zi+i - Zi, z i = Uzi + z i+ i). (2.1) 

*+2 Azi *+2 2 

The resulting shear rates are shown as solid points in Fig. |4j 

In Fig. ^ a) , we observe that the shear rates naturally fall into three sections, marked by the 



horizontal dashed lines in both Fig. [3] and Fig. |4j The uppermost (z > 1) and lowermost (z < 0) 
regions are the boundary layers. When the height of the sample is measured in real units, the width 
of each boundary layer is equivalent to one large particle diameter or two small particle diameters. 
We employ a linear transformation to ensure that z = and z = 1 correspond to the top of the 
lower boundary layer and the bottom of the upper boundary layer, respectively. We limit our 
modeling to z £ [0, 1] since we are interested in the bulk behavior of the system. We also normalize 
the velocity so that u(0) = 1. Since there is no data point at z = 0, we calculate the line containing 
the points (us, z§) and (uq, Zq), which span z = 0, and use the velocity value associated with z = 
on that line to normalize the velocity data (see Fig. 3fta)| ). Note that the velocity of the bottom 



plate sets an overall timescale that is necessary to make a full comparison between predictions of 
the model and the observed segregation in the experiment. However, in this paper we consider 
only a comparison between the theoretical predictions of Case I and Case II, using the experiment 
solely to provide physically realistic shear rate parameters. 

In Fig. |^[bj] we observe that the shear rates can be split into a low-shear region and a high-shear 
region. The division occurs at z c , which we take to be located midway between two adjacent z i 

points: z c = z\q = 0.29. To determine shear rate parameters ko and k\, we average the shear rates 
in each of the two regions. This yields ko = 2.4 ± 0.4 for < z < z c and k% = 0.31 ± 0.05 for 
z c < z < 1, which are both shown as vertical dotted lines in Fig. |^[b)[ For comparison, we can 
use (ko,ki,z c ) to generate the corresponding piecewise linear fit to the measured velocity profile; 
this is shown by the dotted lines in Fig. |3ftb)| These three parameter values are used with the 
constructions of ^3] to generate the solution in Case I shown in Fig. [2^a). 

The measured velocity profile in the region < z < 1 is also well-described by an exponential 



function u(z) = be Z//A + c, as shown in Fig. c [b) A least-squares fit provides model parameters 
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A = 0.22 ± 0.01 and b = 0.82 ± 0.05 and c = 0.14 ± 0.06. The resulting shear rate \u'(z)\ = \eT z l x 
is plotted as a straight line in the semi-logarithmic plot Fig. 4j[a)| and as a curved line in Fig. |^[b)| 



these figures verify that the procedure for determining the exponential shear rate from the measured 
velocity profile also provides a good fit to the shear rates. The parameter values A and b are used 
with the constructions of ^3] to generate the solution in Case II shown in Fig. [2^b). 

To summarize, we have determined parameter values from the experiment for shear rates in 
Case I and Case II. The specific values we use to generate the solutions shown in Figure [2] are: 

Case I: z = 0.5, z c = 0.29, k = 2.4, k x = 0.31. 

(2.2) 

Case II: z = 0.5, A = 0.22, b = 0.82, a = b/X = 3.7. 

3 Initial Boundary Value Problems 



In this section, we derive solutions of the initial boundary value problem (1.2-1.4) in Cases I and 
II. Since the segregation rate parameter s > simply affects the time scale, we first set s = 1, 
and later normalize the time scale by choosing the value for s which provides t* = 1 in Case II. 
We begin with a treatment of characteristics and shocks, focusing on differences from standard 
constructions. 

3.1 Characteristics and Shocks 

Characteristics reduce the construction of continuous solutions of scalar first order PDEs to solving 



ordinary differential equations. For equation (1.2) with s = 1, characteristics are curves z = z(t) 
and if = <p(t) given by 

§ = ^ = -a>(z)f(<p). (3.1) 

Thus, a(z)f((f) is conserved along characteristics: 

a(z)f((p) = constant. (3-2) 

Along characteristics in Case I, in which a(z) is piecewise constant, z(t) is piecewise linear with 
a jump in slope across z = z c and ip is piecewise constant with a jump across z = z c . In Case II, 
the characteristics are smooth curves: since a'(z) ^ 0, the only characteristics which are straight 
lines are those with ip = or <p = 1. All other characteristics are not straight, and moreover, ip is 
not constant along them. This is in agreement with the observation that <p = and ip = 1 are the 
only constant solutions of the PDE in Case II. 

Shock waves satisfy the Rankine-Hugoniot condition, in which the speed of the shock is related 
to the flux across it. Specifically, if the shock is z = j{t), and (p±(t) = ip( r y(t)±,t),a±(t) = a(7(i)±) 
are the one sided limits, then 

dl = a + (t)f(v+(t))-a-W(<P-(t)) ^ (33) 
dt — (P-it) 

This formulation is consistent with the interpretation of the interface z = z c in Case I as a stationary 
shock, for which -yit) = z c , and across which the fluxes balance: 

k f(<p-(t)) = k 1 f(<p+(t)). (3.4) 
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Figure 5: The flux balance condition (3.4) in Case I. 



The flux balance (3.4) is also consistent with the structure (3.2) of characteristics. In Fig. [5] we 
show typical fluxes in Case I with values of cp± representative of the solution of our specific initial 
value problem. In the figure, f((f) has a minimum at ip = (p m . In our case, ip+ < ip m is known, and 
</?_ < p + is then determined from (3.4). However, if </?_ < ip m were given and /cq/(9?-) < kif((p m ), 



then there would be no value of 99+ satisfying (3.4). Consequently, the solution of the initial value 
problem would be rather different, with a shock wave reflected from the interface z = z c . 



3.2 Case I 



In this subsection, we solve the initial boundary value problem (1.2-1.4) in Case I, in which a(z) 



is given by (1.5). For the solution, it is crucial that k$ > k\, which is the physically-meaningful 
relationship for granular shear bands. In addition, since z$ = ^ in the experiment, we also assume 
z* = 1 — zq > z c . In Fig. [2^a), we show the solution with values of ko and k\ calculated from the 



measured velocity profile in §2.1| 

We construct the solution (p(z, t) in several steps, corresponding to the different features in 
Fig. [2^a). For small t > 0, the solution consists of a single rarefaction wave centered at (z,t) 
(zo,0). The rarefaction reaches the material interface 



at a time t = t c , and is transmitted 
through through the interface as a simple wave, in general not centered. The simple wave first 
reaches the boundary z = at a time t = to, and the rarefaction wave first reaches the boundary 
z = 1 at a time t = t\. The simple wave and rarefaction are reflected from the boundaries as shock 
waves z = 7j(i), t > tj, j = 0, 1. The shock z = jo(t) crosses the interface z = z c at a time t = t c , 
and meets the shock z = 71 (t) at a time if. For t > t* , the solution is the piecewise constant 
function 

[l, z < z* 
<p{z,t) = I 

[0 z > z*, 

where z* = 1 — zq, as expected from conservation of the total mass (or volume) of small particles. 
Note that if z* < z c , then the descending shock z = Ji(t) reaches the interface z = z c before the 
rising shock z = 70 (t), and is transmitted through the interface; apart from this difference, the 
solution is the same. 

To simplify some of the construction, and carry explicit calculations as far as possible, we 
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restrict attention to the case /(</?) = <p(<p — !)• Then, the centered rarefaction is given explicitly by 



f = <pi{z,t) 



1 ( z- z 
kit 



+ 1 



z c < z < 1, t > 0. 



(3.5) 



It reaches the boundary z = 1 at time t% = (1 — zo)/ki, since the first characteristic to reach this 
boundary carries <p = 1. The rarefaction is reflected as a shock z = 71 (t) satisfying the jump 
condition (3.3). Since the shock has the centered rarefaction on one side, and i/? = 0on the other, 

71 



it is determined from the initial value problem 

7i zo 



dt 



2t 2t 



h 
2 ' 



1. 



(3.6) 



Thus, 



jl(t) = z - kit + 2V(1 - zojht- 



(3.7) 

Similarly, the characteristic with tp = reaches z = z c at time t c = (zq — z c )/ki. Along the line 
z = z c and (p takes values 



<P+(t) 



1 



t c 



(3. 



The line z = z c in the (z, t)-plane behaves as a stationary shock as far as the weak solution is 
concerned. Consequently, (p jumps from p + to a value <P-{t) = p(z c —, t), t > t c , while keeping the 
flux continuous; the jump condition ( |3.4[ ) is 

k (p-((f- - 1) = kip + (ip + - 1). (3.9) 

Solving this quadratic equation for € [0, |), we find 



¥>-(*) 



(1 



i) 



with = ip+(t) given by (3.8). 

Next, we construct the simple wave that emanates from the line z = z c 
involves a family of straight line characteristics parameterized by r > t c : 

z = kv(2(p — l)(t — t) + z c , t > t. 



(3.10) 



The construction 



(3.11) 



On each characteristic, p = P-(t) is constant. Thus, equations (3.8), (3.10), (3.11) define p(z,t) 
implicitly in the simple wave. 

In order to calculate the shock wave z = 70 (t) that reflects from the boundary z = 0, we need 
to be able to calculate <p(z, t) in the simple wave. Apart from the outermost characteristic 



z = -k (t - t c ) + z c , 



(3.12) 



on which <p = is constant, we find p(z, t) numerically by solving a quartic equation, derived as 
follows. 

The function P-(t) is defined by (3.8), (3.10). We can write the inverse of this function, 
obtaining r = r(ip) : 



r(<p) 



fci 



fc (2 ¥ ?-l) 2 - k + ki 



tc- 



(3.13) 
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Substituting into equation (3.11), we have an equation denning ip as a function of z and t. Let 

t c > 0. Then in the new parameters and variables, 



if) = 2ip - 1, a = 1 - {j- > 0, and /3 



d3Hb, d37T3fc become 



z = koip t 



ft 



VW 2 



+ z c 



a 



(3.14) 



Rearranging and expanding, we find that we have a quartic equation for ip: 

g(ip; z, t) = A^ + Bip 3 + Cip 2 + Dip + E = 0, 
with coefficients depending on z,i given by 

A = kit 2 , B = 2k t(z c - z), C = (z c - zf - kl(at 2 + /3 2 ), D = -2k t(z c - z)a, E = -a(z c - z) 2 . 

(3.15) 

For (z, t) in the simple wave, we seek to solve equation (3.14) for -0 £ (—1, —y/a), corresponding to 
< if < — y/&)- The solution can then be used to find the shock wave z = 70(f)- 



Lemma 1 For (z,t) in the simple wave, g(—l;z,t) > > g(—y/a;z,t), with g(—l;z,t) = only 
on the characteristic (3.12). 

Proof: It is straightforward to check g(—l;z,t) = on the characteristic ( |3.12 ), so we suppose 
that (z, t) lies above that characteristic in the {z, t) plane. First we will show g{— 1; z, t) = A — B + 
C — D + E > 0. Substituting in the values for the coefficients and simplifying, we find 

g(-l; z, t) = kohit 2 - t 2 c ) + 2ht(z - z c ) + ^(z c - z) 2 . 



Then we substitute for z using equation ( |3.11 ) and simplify, concluding that 
g(-l; z, t) = k h (r 2 - t 2 c + 4ipr(t - r) + Aip 2 (t - r)) . 



Along the characteristic (3.11) in the simple wave, we have t > r, and r > t c except along the 
straight characteristic along which <p = 0. Therefore, g{—l;z,t) > 0. 
At i[> = —\J~ol, a similar calculation yields 



g (-- v/«; z, t) = -h (k - h) t 2 c < 0, 
since ko > k\ > 0. This completes the proof. ■ 



Corollary 1 For (z, t) in the simple wave, g(ip; z, t) = has a solution in the interval [—1, —yfa) 
and a positive solution. 



Proof: The Lemma establishes the solution in [—1, —y/a). Since the constant E in (3.14) is neg- 
ative, the product of the four solutions of g(ip; z,t) = is negative. Thus, whether the polynomial 
has all real roots, or two real and two complex conjugate roots, at least one of the roots must be 
positive, since we already have established a negative root. ■ 
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The leading edge of the simple wave is the characteristic (3.12) on which <p 
the boundary z = at time t = to given by to = t c + j^. From the point (z, t) - 
z = 7o(t) emerges from the boundary. Behind the shock is a layer of small particles, with (p 



- 0. It reaches 
(0, to), a shock 



1. 



Consequently, from the Rankine-Hugoniot condition (3.3), the reflected shock satisfies 



dr/o 
dt 



kof(lo,t), 7o(io) = 0, 



(3.16) 



where <p(jo,t) is the value of 92 in the simple wave at the shock. 



Equation (3.16) is solved numerically, since we do not have a closed formula for the simple 



wave. The Corollary shows that <p(p(o{t),t) can be determined by solving equation (3.14). To solve 



equation (3.16), we therefore use the MATLAB function roots in conjunction with the MATLAB 
routine ode45, employing the values (2.2) determined from the experimental data in £2.1 At each 
call of roots, we verify that g(tp;z = 7o(t),t) has two complex roots, thereby checking that we 
have found the only relevant value of 92 in the simple wave. 



As a further check, we compare the coefficients in equation (3.14) with an established criterion 



for the existence of just two real roots. To do so, we place the quartic equation into a normal form 



x 4 + px 2 + qx + r 



0. 



(3.17) 



by dividing (3.14) by the coefficient A, and letting x = ip + tj. The coefficients p,q,r are then 
functions p, q, f of (z, t), in addition to the parameters ko, hi, z c , Zq. Equation (3.17) has coincident 



roots on the swallowtail surface S generated by eliminating x from equation (3.17) and the equation 



Ax s + 2px + q = 0. 



A convenient parametrization of S is obtained by expressing (p, q, r) in terms of p and x : 



-2px 
-qx — px 



4x 3 
.2 



x 



■ px 2 + 3x . 



(3.18) 



(3.19) 



For the parameter values (2.2), 
< z < z c ,t > t c and (q,r)(z c ,t) 



we easily verify that p(z, t) < for (z, t) in the simple wave 
= 0. Moreover, the surface S = {(p, q, f)(z, t) : < z < z c , t > t c } 



lies below the swallowtail S. This region in (p, q, r)-space corresponds to coefficient values for which 



(3.14) has exactly two real roots. In Fig.rol we show the projection of the swallowtail onto the (q, r) 



plane for values of p including the range of p. We superimpose the corresponding projection of S. 



Returning to the structure of the solution of the initial boundary value problem, we observe 
that, although the shock is initially tangent to the i-axis, it then immediately has positive speed, 
since ip(z,t) > for t > to, z > in the simple wave. Consequently, the shock z = 7 o(t) reaches 



z = z c at a finite time t = t c . Since ip± = 1 satisfies the compatibility condition (3.9), the shock 
z = 7o(i) is simply transmitted through the interface but now satisfying the ODE 7 (t) = ^192(70, t), 
with initial condition 70 (t c ) = z c , and ip = (p±(z = 70, t) given by the centered rarefaction wave 



(3.5). Solving the initial value problem, we find an explicit formula for the solution 

7o(*) = zo + ht+ < t(z c - z - hQ, t>t c . (3.20) 
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-0.15 



Figure 6: Projection of the swallowtail and coefficients of g(ip;z,t) in the simple wave in Case I, 



showing the number of real solutions of equation (3.17). The figure on the right is a magnification 
of the dotted rectangle in the left figure. 



To determine the time t* at which shocks 70 and 71 meet, we first note that by mass conservation, 
they meet at the location z = 1 — zq. Then t* can be determined from the equation 71 (t*) = 1 — zq, 
resulting in the expression 

t* = ^ (1 + 2^0(1-20)). (3.21) 

But then 70 (t*) = 1 — zq becomes an equation for t c , with the result that t c is independent of ko, 
and can be calculated explicitly, without resorting to the numerical values of the shock z = 70 (i) 
as it approaches z = z c from below: 

t c = ^ (v^o + v 7 ^) • (3-22) 
In Fig. [2^a), we show i c and t* normalized by the segregation rate constant s = 13.60. In the 



figure, we use the parameter values (2.2). Then (3.21), (3.22) give 

Case I: t c = 0.37; t* = 0.47, 

in agreement with the simulation. 

It is remarkable that these two times are independent of the shear rate ko within the shear band. 
However, this is a simple consequence of conservation of mass. The time t c is the time at which 
enough small particles have dropped below z = z c to form a layer of small particles of depth z c . 
These particles necessarily are transported from z > z c , where their dynamics are independent of 
ko- More precisely, the conservation law cpt + [a{z)f{(p)) z = implies, together with the boundary 
conditions, that cp(z, t) dz is independent of time. But Jq <p(z, 0) dz = 1 — zq, and for t > i c , we 
have ip(z, t) = 1, < z < z c , so that 

1 — zq = / ip(z,t)dz= / tp(z,t)dz+ / ip(z, t) dz = z c + / ip(z,t)dz. 

JO JO J z c J z c 

Consequently, t = t c is the first time for which 

m{t) 

if(z, t) dz = 1 — Zq — Z c , 
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Figure 7: Characteristics (3.1) in Case II. (a) Phase portrait, (b) Projection onto (t, z) plane. 



Curves are numbered by values of ipo in (3.28). 



an equation that does not involve ko, and which can be solved explicitly for t 
f(ip) = <p(tp — 1). Since t c is independent of ko, it follows that t* is as well. 



t r in the case 



3.3 Case II 



First we examine the structure of the solution to (1.2-1.4) for smooth functions a(z) and fluxes 
f{<p), satisfying the following conditions, consistent with the specifications of Case II: 



/"O)>0, /(0) = /(l) = 0, a(z)>0, a'(z)<0. 



(3.23) 



Under conditions (3.23), the invariance (3.2) of a(z)f(ip) along characteristics is easily visualized, 
and gives the phase portrait for the vector field ( |3.1| ), shown in Fig. [7] (using f((f) = tp{tp — 1), and 



parameter values (2.2)). 



From (3.23), f(<p) has a unique minimum, at ip m G (0, 1 
at <p m , and decreasing in p, as shown in the figure. We also have f'(<p) > for ip > ip 



Trajectories of (|3.1|) are horizontal 

m , so that 

(3.1), (3.23) imply that z(t) is increasing there, and decreasing for p <C tpm* 

For fixed zq £ (0, 1), the characteristic curves through (z,t) = (zq,0) form a fan, with <p = 
corresponding to the line on the z axis in the phase portrait of Fig. [7]^a), and similarly ip = 1 
corresponding to the vertical line <p = 1. As remarked earlier, these are the only straight character- 
istics, and the only characteristics on which ip is constant. In between, the characteristics approach 
z = monotonically if <p> < ip m , or have a maximum in z before reaching z = 0, or reach z = 1 
monotonically, before the maximum is reached. This fan of characteristics forms the rarefaction 
wave emanating from (z,t) = (zq,0), and joining ip = to <p = 1. In Fig. W[h), we show the 
characteristics in the (z, t) plane, and remark that the pattern of characteristics is quite different 
from the pattern of contours of ip in the representation in Fig. [2^b) . 

The rarefaction is reflected from the boundaries at z = 0, z = 1, forming a pair of shock waves 
z = ji(t),£ = 0, 1 that eventually meet at a time t = t* , after which the solution consists of the 
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single stationary shock from p = 1 to ip = 0. Since f(p) is smooth and is zero at p = 0, 1, let 
f{p>) = <p(<p — l)g(tp), where g(<p) is smooth and positive. Since <p = 1 behind the shock 70, and 
<p = behind the shock 71, the ODE (3.3) become 



^7o ^7i 

— = 0(70)^(70,^5(^(70,*)), -jr = 0(71) (95(71, *) - 1)5(^(71, *)) (3-24) 

In the general case, the values of p in the rarefaction are known only implicitly, and the shock waves 
can be determined only numerically. Even when these rarefaction values are given by a formula, as 
in the example presented below in connection with the experiment, the ODE may be intractable 
to explicit solution. 



3.3.1 Exponential Shear Rate 



We are able to calculate the explicit solution, except for the reflected shocks, if we tak e a( z) 
a§e~ z l x and f(ip) = tp(<p — 1). In Fig. b[b), we show the solution using parameter values (2.2). 



In the rarefaction wave, tp is defined implicitly by (3.2), which simplifies to the quadratic 
equation 

(3.25) 



<p(l -p) = <p (l - p )e^- z "^\ 



On the curves (3.25), one for each <p £ [0,1], <p and z evolve in time according to (3.1), with 
initial conditions p(0) = po an d z(0) = zq. In particular, the evolution of (p is independent of the 
evolution of z, and in fact, p decays linearly along characteristics: 



dp 
~dt 



a(z) 



p(p - 1) 



a 



(1 - p o )e- Z0/X , <p(0) = p . 



Integrating yields 



p(t) = -a'(z)p {l - p o )e- Z0/x t + p . 



(3.26) 



(3.27) 



Substituting into equation (3.1) gives the evolution of z: 

dz 



dt 

Therefore 



a(z)(2(p-l) = a e- z/x (2Vt + 2^-1), z{0) = z ; V 



a 



<p (l - p o )e- Z0/X . (3.28) 



a 
A 



{Vt 2 + (2p Q - l)t) 



A 



Oo 

A 



p (l - <p )e- z °l x t 2 + 2p> Q t - tj . 



(3.29) 



Eliminating po between (3.27) and (3.29) gives ip = (p(z,t) in the rarefaction wave. Specifically, the 
quadratic equation (3.29) can be solved for p a = <p (z, t). Then, p = p>(z, t) is obtained from (3.27), 
with p = p (z, t). The characteristics in Case II using the experimentally determined parameters 
are shown in Fig. [7| 

Before specifying p(z,t) completely, we consider the leading edges of the wave emanating from 
(z,t) = (zo,0). The leading edge z = z m i n (t) of the rarefaction approaching the lower boundary 
z = 0, carries p = <p Q = and is the particle path of the first small particle to reach the lower 
boundary. Similarly, the leading edge z = z max (t) of the rarefaction, on which p = p a = 1, 
approaches the upper boundary z = 1. Let to, h be the times at which these curves reach z = 
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0, z = 1, respectively, so that z m i n (to) = 0, z max (ti) = 1. These times are easily found from 



equation (3.29) by substituting the relevant values for z and ip - 

A 



k = — I e 
a Q 



z /X 



1 



A 

h = — l r 
a 



l/A _ z /A 



(3.30) 



Now we can solve the quadratic equation (3.29 ) for ip = ip (z, t), using the fact that z = z max {t) 



should give ip a = 1, to select the correct root of the equation: 



^ o(M) = 2aot ( X " 2g2 ° /A + ^ 4e(Z+ZQ)/X + ( 



(3.31) 



The entire rarefaction fan is now characterized using equation (3.27) with tp given by (3.31): 

(3.32) 



<p(z,t) 



A 



t<p (z, t)(l - <p (z, t))e- Zo/X + <p {z t t). 



Next, we formulate an ODE for the reflected shocks z = 7o(i)> z = 7i(t), originating from 
(z,t) = (0,to), and from (z,t) = (l,ti) respectively. In order to track the shocks using the 



Rankine-Hugoniot condition, we use the expression for ip = (p(z,t) given in equation (3.32) in the 
region between the shocks and between the outermost characteristics z = z m i n (t) and z = z max (t), 
that is, the region in which we find the rarefaction fan. 



As in Case I, the Rankine-Hugoniot condition for a shock curve z = j(t) for the PDE (1.2) is 

(3.33) 



given by (3.3). Therefore, for 7 = 7o(i), where <p jumps from 93 = 1 to 99 = <p(jo,t), we have 

dr/o 



dt 



a(7oM7o, t), t > t , 7o(*o) = 0. 



Similarly, for 7 = 71 (t), where (p jumps from ip = 95(71, t) to ip = 0, we have 



dt 



a(7i)(¥>(7i,<)-l), * > h, 7i(*l) = 1- 



(3.34) 



We solve the equations (3.33) and (3.34) using the Matlab ODE solver ode45. The time t* at 



which 70 and 71 meet is calculated numerically. The entire solution is shown in Fig. |2jb), where 
the time scale has been normalized by the calculated t* = 13.60; this is equivalent to setting the 
segregation rate s = t* . 



4 Discussion 

In both cases, the structure of the solutions in ^3] captures the segregation process observed in 
the experiment, but there are significant differences between the two solutions, and between the 
solutions and the experiment. 



Having computed solutions in each of Case I and Case II using the parameter values (2.2) 
derived from experimental data, and having noted the differences in construction, we are left with 
the observation that the time to segregation in the two cases is markedly different, with the ratio 
of final times t* being approximately 0.47. Because the exponential shear rate of Case II is a better 
fit to the data, it is natural to regard it as the better model. 
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There are additional reasons to favor Case II over Case I. It would be tempting to blame the 
values of the shear rates k® and k\ for the lack of agreement in t* , but t* is in fact independent of ho- 
Therefore, the problem has to lie in the upper part (z > z c ) of the sample. Indeed, the value of t\ in 
Case I is significantly smaller than in Case II, and moreover, the shock 71 in Case I is steeper over 
most of its path than the corresponding shock in Case II. Both of these effects would be countered 
by decreasing the shear rate k\ so that small particles approach the interface z = z c more slowly. 
However, the answer may be more subtle. In Case I, k\ overestimates the experimental shear rate 
(and that of Case II) near z = 1, thereby promoting segregation there. This appears to be a more 
significant effect than the retarding of segregation closer to z = z c due to k\ underestimating the 
shear rate there. It is also significant that, as indicated by the structure of characteristics in Case II 
(Fig. [7]), the upper part of the domain strongly influences the segregation in the lower part, so that 
in Case II, the entire domain is involved in determining the time to segregation. Consequently, the 
larger segregation rate near z = has significance in Case II but none at all in Case I. In summary, 
while the broad structure of the solutions (rarefaction wave and reflected shocks) is captured in 
both Case I and Case II, any quantitative comparison to experiments must utilize the more refined 
smooth fit to the shear rates of Case II. 

The comparison to experimental results faces a number of challenges, discussed in greater detail 
in our companion paper [9]. Here we note: 



1. The model does not account for the rapid opening of void spaces due to Reynolds dilatency [H] 
during the initial transient dynamics. During this process, particles exhibit a more disorderly 
motion and the resulting measured velocity profile is inconsistent with those measured at 
later times. This short-time behavior undoubtedly accelerates the initial mixing of small and 
large particles prior to the establishment of the conditions described by the model. Therefore, 
the model only applies after this initial transient period, by which time the configuration of 



the particles is no longer given by (1.3) and is instead a somewhat mixed state. 



2. Full segregation is never achieved in the experiment, making it impossible to identify a final 
time t* to compare with the models. To measure the progress of segregation in the experi- 
ment, we record the height of the top plate and relate the total volume of the sample to the 
degree of mixing/segregation. Because we observe that the volume approaches its final value 
only exponentially in time, we conjecture that isolated large particles remain trapped at the 
bottom and migrate upward on a slower timescale. This is a discrete effect not captured by 
a continuum model. 

3. The model does not account for three-dimensional motion of the particles. While the as- 
sumption that the concentration of large and small particles depends only on depth is clearly 
unrealistic, (f(z, t) might be considered reasonable as an average across horizontal layers of 
particles. The side walls of the experiment undoubtedly influence the dynamics, since they 
introduce lateral shear. It would be possible to include this effect in a multidimensional 
model, but verifying such a model with experiments would be hampered by the difficulty of 
tracking particles within the bulk. 

Finally, we note that the piecewise constant case (Case I) has been studied both theoretically 
and numerically in various contexts, for example in flow in porous media composed of layers of 
different material such as sand and clay, and in connection with sedimentation [TJ El [TTJ, [18] . In 
these studies, a variety of techniques are introduced for studying entropy conditions and special 
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solutions, as well as analysis of existence and uniqueness questions. Our application to segregation 
solves a special initial boundary value problem, but it reveals the unforeseen consequence that the 
material interface at z = z c removes the influence of the lower portion z < z c of the domain on the 
overall time scale of the evolution. 
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